#########
#Replication Code: "Policy Feedback, Energy Equity, and Climate Justice: Can Existing Policies Improve Solar Access for Low- and Moderate Income Communities in the U.S.?" 2024 PS: Political Science & Politics
#Aparajita Datta, University of Houston 
#adatta3@central.uh.edu"
#Last Update: July 1, 2024
#######
  

#Require packages, define functions, and set working directory. Note: Include name of working directory in line 11. 
rm(list=ls())
setwd("include name of directory")
require(maxLik)
library(AER)
library(foreign)
library(ggplot2)
library(MASS)
library(gmodels)
library(Hmisc)
library(MASS)
library(ordinal)
library(reshape)
library(stats)
library(ggplot2)
library(KMsurv)
library(nlme)
library(simPH)
library(survival)
library(car)
library(flexsurv)
library(readxl)
library(plyr)
library(stargazer)
library(plm)

#Read data from excel file. Note: Include the name of the file in line 34. 
my_data <- read_excel("include file name")

#Descriptive analysis of the variables. Available in Table S1 of Supplementary Material. 
summary(my_data)

#Recode party data for state legislatures and governors as numeric 
my_data$gov <- revalue(my_data$Governor, c("Rep"= "1", "Ind"="2", "Dem"=3))
my_data$LControl <- revalue(my_data$LegislativeControl, c("N/A"= "0", "Rep"="1", "Split"="2", "Dem"="3"))
my_data$gov_num <- as.numeric(my_data$gov)
my_data$lc_num <- as.numeric(my_data$LControl)

#Event History Analysis using Cox Proportional Hazards Models
#Model 1 - Efficiency and Electricity Sector Portfolio Standards only 
model1 <- coxph(Surv(start,end,Event) ~Efficiency + ElectricitySector,data= my_data)
summary(model1)

#Model 2 - Share of Neighboring States only
model2 <- coxph(Surv(start,end,Event) ~ NeighboringStatesEfficiency + NeighboringStatesElectricitySector,data= my_data)
summary(model2)

#Model 3 - Share of Neighboring States and Federal Spending only 
model3 <- coxph(Surv(start,end,Event) ~ Spending + NeighboringStatesEfficiency + NeighboringStatesElectricitySector,data= my_data)
summary(model3)

#Model 4 - Full Model 
model4 <- coxph(Surv(start,end,Event) ~Efficiency + ElectricitySector + Sunshine + SolarShare + Bill+ InterestGroups + Poverty + Spending + gov_num+lc_num+ NeighboringStatesEfficiency + NeighboringStatesElectricitySector,data= my_data)
summary(model4)

#Present results as in Table 1
require(stargazer)
stargazer(model1,model2, model3, model4, type = "text", dep.var.labels = c("Time to LMI adoption"), covariate.labels = c("Efficiency Policy", "Electricity Sector Portfolio Standards", "Solar Potential", "Share of Solar in Energy Mix", "Average Annual Electricity Bill", "Renewable Energy Interest Groups", "Share of Population below FPL", "Federal Spending on Renewable Energy", "Governor's Party", "Legislative Control", "Neighboring states- Efficiency", "Neighboring states- Electricity Portfolio Standards"), out = "fit1.html")

#Schoenfeld residuals diagnostic test. Available in Table S2 of Supplementary Material.
z <- cox.zph(model4)
par(mfrow=c(2, 2))
plot(cox.zph(model4))
z


#Figure 3 Survival function for Model 4- Full model
pdf(file="survivemod.pdf",heigh=6,width=6)
plot(survfit(model4), ylim=c(0.5, 1), xlab="Years",
     ylab="Proportion of states - not adopted LMI incentives")
dev.off()

#Supplementary Material - Battery of dummies - Table S3
modelS3 <- coxph(Surv(start,end,Event) ~factor(Efficiency) + factor(ElectricitySector) + Sunshine + SolarShare + Bill+ InterestGroups + Poverty + Spending + NeighboringStatesEfficiency + NeighboringStatesElectricitySector,data= my_data)
summary(modelS3)

#Supplementary Material - Dynamic Logistic Regression - Table S4
modelS4 <- plm(Event ~ Efficiency + ElectricitySector + Sunshine + SolarShare + Bill+ InterestGroups + Poverty + Spending + gov_num+lc_num+ NeighboringStatesEfficiency + NeighboringStatesElectricitySector,data= my_data, index="TimetoLMI", model ="random", random.method ="walhus")
summary(modelS4)

#Replication code ends here
